Optical Forces Between Coupled Plasmonic Nano-particles near Metal Surfaces and 

Negative Index Material Waveguides 
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We present a study of light-induced forces between two coupled plasmonic nano-particles above 
various slab geometries including a metallic half-space and a 280-nm thick negative index material 
(NIM) slab waveguide. We investigate optical forces by non-perturbatively calculating the scattered 
electric field via a Green function technique which includes the particle interactions to all orders. 
For excitation frequencies near the surface plasmon polariton and slow-light waveguide modes of the 
metal and NIM, respectively, we find rich light-induced forces and significant dynamical back-action 
effects. Optical quenching is found to be important in both metal and NIM planar geometries, which 
reduces the spatial range of the achievable inter-particle forces. However, reducing the loss in the 
NIM allows radiation to propagate through the slow-light modes more efficiently, thus causing the 
light-induced forces to be more pronounced between the two particles. To highlight the underlying 
mechanisms by which the particles couple, we connect our Green function calculations to various 
familiar quantities in quantum optics. 

PACS numbers: 42.50.-p, 78.67.Bf, 73.20.Mf, 78.67.Pt 



I. INTRODUCTION 

Since the proposal of using intense laser light to trap 
atoms 1 and particles 2 by Ashkin, scientists have achieved 
a remarkable ability to manipulate matter via optical 
forces''. This has lead to a plethora of light- matter force 
manipulation techniques, from Gaussian beam optical 
traps , that are now routinely used in laboratories', to 
near-field nanometric tweezers'' which have motivated an 
entire field involving plasmonics to enhance and manip- 
ulate optical forces '~ n . The use of plasmonic structures 
allow strong evanescent field enhancements near the plas- 
mon resonances of the system, and can efficiently trap 
particles with a size smaller than the wavelength of illu- 
mination. Examples include a patterned substrate with 
metallic particles 10 ' 11 , metallic near-field tips' 1 , or close to 
nano- apertures in thin metallic films 12 . In the latter case, 
Juan et. al. 12 recently examined the trapping of 50 nm 
and 100 nm polystyrene particles near a nano-aperture in 
a thin sheet of gold and illuminated with light very near 
to the transmission cutoff wavelength of the aperture. It 
was observed that as the particles were trapped in this 
geometry, the scattering induced by the particle worked 
to enhance the trapping, causing a self-induced "back- 
action" . This back-action illustrates that when consider- 
ing the optical forces on nano-particles (NPs) near reso- 
nances, the NPs can interact non-perturbatively with the 
system and thus any models or theoretical descriptions 
must include NP coupling in a self-consistent way. 

Similar to field-enhancements using plasmonic struc- 
tures, it is also possible to obtain radiation enhancements 
using metamaterial structures which can be engineered 
with a negative index of refraction (e < 0, /x < 0). Vese- 
lago introduced the concept of a negative refraction in 
1968 1! , and, in 2000, Smith et. al. experimentally real- 
ized such a material. This breakthrough initiated the 
field of "metamatcrials" in which patterned materials 



with unit cells smaller than the wavelength of illumina- 
tion were created to tune the effective material responses. 
Soon after, predictions were made of super- lensing 1 5 , 
cloaking devices 1 ' 3 ' 1 ' and "left-handed" waveguides 18 , all 
prompted by the ability to create negative index materi- 
als (NIMs). 

Negative index materials also posses interesting quan- 
tum optics and quantum electrodynamics (QED) prop- 
erties. For instance, Kastel and Fleischhauer 1 " exam- 
ined an idealized (non-absorbing) negative index mate- 



E 








FIG. 1. (Color online) Schematic of the geometry. We con- 
sider the optical forces on two NPs with various resonance 
frequencies and radii of 0.015A in air above a planar geom- 
etry. We study both an infinite half-space of silver, and a 
280 nm metamaterial slab which supports slow light modes 
(see text for details). The NPs are illuminated with a plane 
wave perpendicular to the interface, with the electric field 
polarization along the axis between the particles. 
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rial placed on top of a mirror with an atom above the en- 
tire structure. They found that depending on the height 
of the atom it was possible to obtain complete suppres- 
sion of spontaneous emission. Additionally, they exam- 
ined atoms separated by a perfect negative index slab 
(n = — 1) and showed that it was possible to obtain per- 
fect subradiance and supcrradiance. Yao et. at 2 " exam- 
ined the enhancement of the spontaneous emission rate 
(Purcell Factor-' 1 ) above a NIM slab which supports a 
collection of slow-light modes in the region of negative 
index 1 *- 22 ; related Purcell factor enhancements in NIM 
waveguides have been found by Xu et al. and Li et 
al. 2i . In the slow light region, very large field enhance- 
ments can be obtained similar to the enhancements found 
in plasmonic structures. However, there were two impor- 
tant differences between the plasmonic and NIM systems: 
(i) the effects of quenching through non-radiative energy 
transfer, and (ii) the role of the quasi-static approxima- 
tion, i.e., the neglect of dynamic retardation effects. In 
plasmonic systems, it is possible to obtain strong Pur- 
cell factor enhancements however most of the emission 
is absorbed by the metal manifesting in non-radiative 
quenching 2 '. This optical quenching is a result of the 
intrinsic material loss which, in principle, can be tuned 
in metamaterial systems. For metals, the quasi-static 
approximation is also typically used 2 ' 1 ' 2 ' for small ob- 
jects approaching the surface (kh < 1, where k is the 
wavevector in the background medium containing the ob- 
ject and h is the height above the surface), but for NIMs, 
such an approximation does not necessarily hold even for 
kz ~ 0.03 20 . 

Another potential advantage of NIMs, is that it is 
possible to obtain large Purcell factor enhancements in 
the optical region of spectrum, which also translates to 
larger spatial distances. Recent experimental results 
have shown the possibilities of engineering the sponta- 
neous emission rate using mctamaterials composed of 
silver nano-wires embedded in PMMA. The nano-wire 
structure was shown to exhibit an in-plane anisotropic 
hyperbolic dispersion curve which have been shown to 
exhibit a negative refractive index 2 ' 1 . By exciting dye 
molecules embedded on top of the metamaterial the de- 
cay rate of the dye molecules was shown to be enhanced 
by a factor of 6 at A = 800 nm compared to dye deposited 
on silver and gold films. This was even greater than ex- 
pected when compared to a semi-classical model which 
predicted an enhancement of 1.8 30 . 

In this work, we theoretically investigate the optical 
forces on NPs above metallic and NIM planar geome- 
tries by self-consistcntly including the interaction of the 
NPs with the scattered electromagnetic field. We ex- 
ploit a photon Green function technique that can be 
used to introduce multiple particles within the dipole 
approximation ' 1 or, if required, by using the full cou- 
pled dipole method 52-34 , where the latter discretizes the 
NPs as small polarizable subunits. Green function and 
coupled dipole approaches have been very successful in 
examining light scattering 35 and optical forces' 51 '' 54 in di- 



electric particles located in evanescent fields above glass, 
and for computing optical forces between metallic par- 
ticles in free space 36 and above glass 21 '. Green function 
approaches have also been used to successfully model op- 
tical trapping using near field optics ' . In addition to 
computing the light-induced forces, we also analyze the 
properties of the Green function, both with and without 
the particles, which contains all the key electromagnetic 
interactions of the system, including coupling to surface 
plasmon polaritons (SPP) of the metallic half-space, the 
localized surface plasmons (LSP) of the particles, and 
evanescent coupling to slow light modes (SLMs) of the 
NIM slab. To help clarify the underlying physics, we 
also make direct comparisons with various well known 
concepts in quantum optics, such as the Purcell Factor, 
the Lamb shift, and real and virtual photon exchange; all 
of these effects are relevant for understanding the ensuing 
coupling dynamics between the particles. 

Our paper is organized as follow. In Sec. II we de- 
scribe the theory of light induced-scattering from NPs 
close to multi-layered surfaces, including, II A - the self- 
consistent calculation of the Green function and the elec- 
tric field, and II B - the calculation of the force from the 
total electric field. We exemplify our theoretical results 
for silver half-space in Sec. Ill, and for the NIM slab in 
Sec. IV. In Sec. V, we conclude. 



II. THEORY 

A. Green Function Calculation 

For a spherical particle with radius a, and in the limit 
where fc^a << 1, the "bare" (i.e., no radiative coupling) 
polarizability is given by the Clausius-Mossotti relation, 



a (w) = 47T£ B a 



3 £ M 



e (w) + 2ei 



(1) 



where s(uj) is the particle dielectric constant (relative 
electric permittivity) imbedded in a homogeneous mate- 
rial with a background dielectric constant, sb, assumed 
to be real. Here ks = oj^/eb/ 1 c is the wavevector in the 
background material. It was shown by Draine ''' that in 
order to satisfy the optical theorem, this polarizability 
must be corrected to include the homogeneous-medium 
contribution to radiative reaction: 



a (lu) = 



a (lu) 
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where the self-induction term, M(lu), can be calculated 
exactly for a spherical particle 38 ' 39 , and for fcsa << 1 can 
be approximated as M = 2i(afc_e) 3 /9. Alternatively, this 
term comes naturally from the homogeneous contribution 
of the Green function. 

Our system is initially characterized by an initial 
electric field E*- ) (r;u>), and an initial Green function 
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G w (r, r';w), in the absence of any particles. The field, 
E(o) (r;w), can be of any form we wish (e.g., Gaussian 
beam, plane wave including reflections from surface), and 
we define it as the field prior to adding any scatterers. 
We subsequently introduce TV particles into the system 
where the ith particle is at position r,, with polarizability 
on (oj). The total electric field - excitation field plus par- 
ticle scattered field - can be calculated self-consistently 
from 



EW (rjw) = E<°> (rjw) 



N 



5>i HG (0) (r,r i;w )-EW (r^w) 



(3) 



Similarly, the Green function of the system after adding 
N particles can be calculated via the Dyson equation 40 , 



G W (r,r' W ) = G (0) (r,r» 



N 



+ ^a 4 (w)G (r,r i; cj)-G (r.;,r';w), 

1=1 

(4) 

which is now the total Green function of the medium, 
including the response of the NPs. These equations form 
the basis of the coupled dipole method 52 , and, impor- 
tantly, they apply to any general inhomogencous and 
lossy media. 

To help clarify the underlying physics we will fur- 
ther assume particles with a size much smaller than the 
wavelength, and consider each NP within the dipole ap- 
proximation; however it should be noted that at very 
short inter-particle distances this approximation eventu- 
ally breaks down 41 . Thus we will restrict the distances to 
regimes where the dipole approximation is expected to be 
a good approximation; in this way, we include the particle 
dipoles exactly, while essentially dealing with spatially- 
averaged particle quantities. Using the Dyson equation, 
we can also rewrite the right hand side of Eq. (3) to be 
given only in terms of the excitation field'*, E^ ) (r;w). 
One has 



EW (r;w) = E<°> (r;w) 



N 



+ £)o i (w)G W (r,r i; a,)-E(°) (r 4 ;w) , 



(5) 
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where G includes the particle(s) response. The 
Green function, G, can also be separated into homo- 
geneous (direct), Ghom, or scattered (indirect), G sca tt, 
contributions. Since Re[Ghom(r, r' — s- r)] diverges, in 
what follows below wc will consider the non-divergent 
Re [G sca tt (r, r')] when r = r', as this is the only rele- 
vant photonic contribution. While Re[Ghom(r, r)] can 
give a very small vacuum Lamb shift, this effect can be 
already included by simply redefining the resonance fre- 
quency of the particles. In addition, since a{u>) includes 
the effect of Im[Ghom(r, r)] through the self-induction 



term, then wc need only consider G scatt (r,r), i.e., when 
r = r'. For r ^ r', wc consider the full G(r,r) = 
G sca tt (r,r') + G hom (r, r'). Frequently, it is possible to 

(TV) 

solve Eqs. (3-5), perturbatively, by considering all G 
and EW on the RHS to be the unperturbed quantities, 

i.e., approximated by G^ and E'°). When the system 
constituents are far separated this can hold; however, as 
the particles come closer together and nearer the planar 
surface, we will show that dynamical coupling become 
important and cannot be neglected. 

To examine the mechanisms of light-induced forces, it 
is useful to consider the Green function of the planar 
surface with and without the inclusion of the particle(s). 
We define the complex local density of states (LDOS) as 



G mm (r,r;w) 
Im[G*°£(r,r;w) 



(6) 



where Gmm 



r,r;w) is given by Eq. (4). Although the 



LDOS only depends on Im[pm '}, we introduce p m as a 
complex quantity for ease of notation. For example, the 
imaginary part of the Green function in a homogeneous 
lossless material at r = r' is given by 



Im[G£r(r,r;< 



IEb 



6ttc 3 



(7) 



which is related to the homogeneous-medium LDOS. 
When the homogeneous material contains loss, then both 
the real and imaginary part of Ghom formally diverge as 
r > r' instead of just the real part of Ghom (i-e., as in 
the case of a lossless material). Consistent with our dis- 
cussions and notation above, when pffl (r; u) = 0, this 
means that the LDOS at that point is equal to the homo- 
geneous density of states, as the scattered contribution 
is zero; thus the total = p m + 1, and we will focus on 
p m . Since we will consider the force interaction between 
two particles, it is also useful to introduce a complex 
non-local density of states (NLDOS), 



d'(r,r» 
Im[G^ m (r',r'; W )]' 



(8) 



which describes light propagation between the two space 
points r and r' 7^ r. 

Many of these quantities are useful also for connecting 
to the quantum optical properties of the optical NPs 42 ' 43 . 
This is a key strength of the Green function approach 
over brute-force numerical electromagnetic techniques 
such as, e.g., FDTD (finite-difference time domain) 44 . In 
the case of the complex LDOS, the real part of Eq. (6) 
can describes frequency shifts (Lamb shifts) of an emitter 
caused by the environment, whereas the imaginary part 
describes the material-dependent spontaneous emission. 
For the photonic Lamb shift, one has 



d ' [G sca tt (r, r ;^)] • d 
he 



(9) 
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for an emitter at position r. For calculations of light- 
induced optical forces, the real part of the LDOS can 
therefore manifest itself as resonance shifts of the local 
electric field. 

The variation of the imaginary part of the LDOS 
causes gradient forces on the particle as it moves through 
the field. In terms of non-local QED interactions between 
two particles, the real part of the NLDOS describes vir- 
tual (instantaneous) photon exchange between two points 
or emitters and the imaginary part describes real (dy- 
namic) photon exchange. Virtual photon exchange man- 
ifests itself in well known processes like Forster coupling, 
which, for a homogeneous medium, exhibit an R~ 3 scal- 
ing with inter-particle separation (R) in free space, and 
real photon exchange manifests itself in dipole-dipole 
coupling which exhibits an dependence 1 '. The effect 
of the real part of the NLDOS on classical optical forces 
manifests itself in the scattered contribution to the force 
and the imaginary part would contribute in a manner 
similar to radiation pressure. For a non-homogeneous 
medium, we stress that the photon coupling mechanisms 
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arc considerably more complicated than simple Forster 
coupling. In addition, we fully include dynamical retar- 
dation effects through the frequency-dependence of the 
response functions. 

The above prescriptions are relatively straightforward 
provided one knows the bare Green functions without 
any particles. For the calculation of the initial planar 
Green function, we use a well established multilayer tech- 
nique 11 ', which is outlined in Appendix A and further nu- 
merical details arc detailed by Paulus et. al. A1 . Although 
we specialize our study for two NPs, the generalization 
to any arbitrary number of NPs is straightforward with 
no change in theoretical formalism. 



B. Light-Induced Forces 

Using the electric field E^* 1 (r, ;w) at the NP position 
r-j, the time-averaged total force on a particular NP from 
a time-harmonic electromagnetic wave is (uj is implicit), 



(F (r0> = § J t (r0 + ^ (r ' i} )*) ' V (r0 + (r * 

+ (aiEW ( ri ) + a* (e^ (r,)) *) x (b W (r 8 ) + (b W (r 8 ) 



(10) 



where £o is the permittivity of free space, B is the mag- 
netic field and T is the period of the time-harmonic ra- 
diation. Upon carrying out the integration, and using 
B = l/iw V x E, and E = -iuiB, 



(ft) = 



>{. 



e jkl e, 



Ei N) d k [E. 



(E^y\ }, 



(ii) 



where e^ kl is the Levi-Civita tensor. Using the relation 



p.jkl _ 
fc '-Imn 
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$in$kmi "we obtain the desired force 
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(N 



(r. t )v(£f° (r.) 



j=x,y,z 



(12) 

This light-induced force describes the force on the parti- 
cle due to the electric field and its interaction with the 
planar geometry as well as the scattering due to the other 
NPs in the system. This is different from the usual gra- 
dient force given by F = ^chqW |E| 2 , which is only ap- 
plicable when a is real and the phase of the field varies 
slowly in space. 



III. SILVER HALF-SPACE 

To calculate the optical forces on the NPs, we consider 
the geometry shown in Fig. 1, with two NPs in air above 



a planar structure. For silver, we consider a half-space 
geometry or an optically thick slab with a permittivity 
given via the Drude model, 



s(uj) 



" j" 



(13) 



Here e r = 6 is the permittivity as to — > oo, ui pe =9.87 eV 
is the electric plasma frequency, and 7 = 51 meV is the 
damping rate due to collisions 49 ' 00 . For a metallic half- 
space, the characteristic surface plasmon polariton (SPP) 
frequency is given by Rc [e (ui ~ u>spp)] = —£b- Below 
this frequency, SPPs are confined to the interface and can 
only be coupled to by breaking the symmetry of the sys- 
tem (e.g., via a grating coupler). For a silver/air interface 
the SPP is located at ujspp = 3.73 eV (Xspp = 332 nm). 
Surface plasmon polaritons are well known for their abil- 
ity to enhance the electric field in their vicinity, but these 
enhancements are also associated with high losses mean- 
ing that coupling between objects or the far field can be 
suppressed / quenched. 

For a spherical NP, the localized surface plasmon res- 
onance is at Re [e (uj)] = —2eb [see Eq. (1)] which occurs 
at wlsp = 3.49 eV (Xlsp = 355 nm) using the param- 
eters for silver and air. To investigate coupling between 
the NPs and the resonances of the metallic half-space, 
we will fix e r and 7 to the parameters for bulk silver de- 
spite deviations from bulk-like behavior for small NPs ' 1 . 
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FIG. 2. (Color online) Particle and silver half-space response 
functions. For both the silver half-space and the silver NP 
permittivity we use the Drude model [Eq. (13)], with 7 = 
51 meV and e r = 6. For the silver half-space, the plasma 
frequency is given by ui pe = 9.87 eV. To tune the LSP of the 
NP to the SPP we set the plasma frequency at uj pe — 10.56 
eV; and to tune the LSP of the NP to be off-resonant with the 
SPP, we set the plasma frequency at ui pe = 7.47 eV. (a) Bare 
polarizability [Eq. (1)] of a silver NP with radius a = 5 nm (= 
0.015 Xspp) with the LSP resonance tuned to the SPP of the 
silver half-space, lospp = ulsp = 3.73 eV (solid lines), and 
with the LSP resonance tuned to be off resonant with the SPP 
of the silver half-space, ljlsp = 2.63 eV (dashed lines). Gray- 
light curve indicates real parts and red-dark curve indicates 
imaginary parts, (b) Refractive index of the silver half-space, 
gray-light curve indicates real part (scaled by a factor of 50) 
and red-dark curve indicates imaginary part. 



For small NPs, the localized surface plasmon of the NP 
becomes strongly dependent on size 52 and shape 53 and 
can be further detuned by adding dielectric or metallic 
coatings 54 . We shall use this as motivation to tune the 
LSP of our NPs to be either on resonance with the SPP 
{ulsp = ^>spp — 3.73 eV, uipe = 10.56 eV), or off res- 
onance with the SPP (u>lsp = 2.63 eV, u> pe = 7.47cV). 
The bare polarizability and metal half-space permittivity 
response functions are shown in Fig. 2. 

In the following, we consider the geometry shown in 
Fig. 1, where NPs with tunable LSPs are located above 
a planar structure. We first use particles with a radius 
of0.015Aspp (= 5 nm) so as to have a reasonable sized 
particle for which the dipole approximation will apply. 
We vary the height, h, of the particles above the half- 
space (measured from the center of the particles) while 
keeping both particles at the same height for simplicity. 
Also, we vary the separation, s, between the particles 
(along the y direction and measured from their centers) 
and the frequency of light illumination. To illuminate the 
particles, we use a homogeneous excitation field, which is 
a solution to the scattering problem (incident light plus 
scattered light) without any NPs; we choose the polariza- 



tion to be along the direction of the particles (y direction) 
to maximize the effect of the coupling between the parti- 
cles. To excite surface plasmons, which only exist for TM 
polarization, the symmetry of the system must be bro- 
ken to allow coupling between the surface plasmon and 
the particles, so only the scattered field from the NPs 
can excite the SPP. The incident intensity is 1 W//xm 2 , 
which we choose only as a convenient reference - the force 
scales linearly with the incident intensity as can be seen 
from Eq. (12). For comparison, we note that the earth's 
gravitational force on a 5 nm silver particle is 54 fN; for 
the intensities considered below, the gravitational force 
is smaller by a factor of 10 -3 , and thus can be safely 
neglected. 

In Fig. 3 (a) we examine the y component of the LDOS 
of the half-space, which will dominate the forces for the 
particular illumination scheme that we have selected. We 
vary the height [r = (0,0, h)] keeping in mind that we 
cannot approach closer to the structure than our particle 
radius which is indicated by the gray shaded region. We 
consider lo ~ lospp — ^lsp arL d add the first particle at 
ri = r and the second particle at r2 = (0, 0.05 Xssp, h). 
When there are no particles in the system, we can see that 
the imaginary part of the LDOS (red-dark dashed curve) 
diverges, which would lead to an infinite LDOS for an ex- 
cited emitter; although this effect may seem surprising, 
such a divergence always happens above a lossy struc- 
ture 20,55 [Eq. (A10)]. However, the inclusion of the par- 
ticle where we calculate the LDOS (red-dark chain curve) 
acts to renormalize the LDOS for distances < 0.08 A. The 
addition of a second particle at the same height as the 
first, but separated by 0.05 A (center to center), further 
renormalizes the LDOS (red-dark solid curve), though it 
becomes apparent that the effect of the silver half-space 
becomes negligible for heights greater than about 0.06 A 
when there are two particles. 

For the real part of the y component of the LDOS 
(gray-light curves), and with no particle in the system 
(dashed), we see that there is a minimum as the par- 
ticle approaches the half-space which would give a blue 
shift for an emitter placed close to the surface. Includ- 
ing a particle at this location (gray-light chain curve) 
reduces the blue shift and including the second particle 
(gray solid curve) causes a change in sign which means 
an emitter would be shifted to the red. For both the real 
and imaginary part of the LDOS these shifts are only 
seen by going beyond the perturbative limit and solving 
Eq. (4) exactly. These results emphasize the pronounced 
back-action effects that occur in describing the electro- 
magnetic properties of the medium. 

Figure 3 (b) shows the yy component of the NLDOS 
in a similar manner to the LDOS described above, with 
lo = lospp = ulsp, r' = ri = (0,0, h), and we consider 
r = r 2 = (0, 0.05 Xspp, h). With no particles in the sys- 
tem, the imaginary part of the NLDOS (red-dark dashed 
curve) becomes large but remains finite as the surface is 
approached, implying that photons are very easily trans- 
ferred from r' to r [off scale). However once a particle is 
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FIG. 3. (Color online) (a) Complex LDOS, p { y N) (r'.w), and 
(b) NLDOS, (r, r',w), for a silver half-space as a func- 
tion of height at ui = uispp- For the LDOS, r = (0,0, ft), 
and for the NLDOS, r' = (0,0, ft), r = (0,0.05 A, ft). Gray- 
light curves represent the real part and red-dark curves rep- 
resent the imaginary part. The dashed line has no parti- 
cles [py (r;w), pyy (r, r';cj)], the chain line has one parti- 
cle located at the position the LDOS/NLDOS is calculated 

[p^ (r;a;), py 1 ^ (r , r' '; u>) and ri = (0,0, ft)], and the solid line 
has two particles - one where the LDOS/NLDOS is calcu- 
lated aand the other at the same height as the first but sep- 
arated by 0.05 A in the y direction [pj, 2 ' (r;w), p y 2 y (r, r';w), 
ri = (0,0, ft), and r-2 = (0, 0.05 A, ft)]. The gray shaded re- 
gion indicates region where particles would overlap the planar 
structure. 



added (red-dark chain curve), the NLDOS reduces dras- 
tically as the particle breaks the symmetry of the sys- 
tem allowing quenching to occur. Interestingly, the addi- 
tion of the second particle (red-dark solid curve) further 
breaks the symmetry of the system and allows light to 
couple to more channels in the planar structure, which 
further enhances the quenching effect and reduces the 
NLDOS. This dramatic reduction is a result of the non- 
perturbative coupling between the particles and the half- 
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FIG. 4. (Color online) lm[p y N ^ (r, to)] for a silver half-space as 
a function of frequency at the position r = (0, 0, 0.05 Xspp)- 
The dashed line has no particles, Im[py ' (r; uj)], the chain line 
has one particle, Im[py (r; to)], and the solid line has two par- 
ticles, Im[py 2 ' (r; w)] at locations ri = (0, 0, 0.05Aspp), and 
r 2 = (0,0.05Aspp,0.05Asfp). (a) The particle LSPs are red 
detuned by — 1.1 eV compared to the SPP frequency 
(see Fig. 2). (b) The particle LSPs are at the SPP frequency. 
Arrows indicate the particular scenarios we are examining in 
force graphs: 'i' - Fig. 5(a), 'ii' - Fig. 5(b), 'iii' - Fig. 5(c) 
and 'iv' - Fig. 5(d). 



space which is theoretically described through the self- 
consistent solution of Eq. (4). For heights greater than 
0.15 A, light propagates purely via virtual photon prop- 
agation as the real part (gray-light curve) approaches a 
finite value (the homogeneous Green function) for both 
zero (dashed), and one (chain) particle. For two particles 
this happens even closer to the surface (ps 0.1 A) due to 
the additional scattering events. Real photon propaga- 
tion occurs when the half-space begins to interact with 
the system and multiple paths are possible for a photon 
to reach r from r'. 

The quasi-static approximation is often invoked for 
particles very close to a surface or to each other 5 ', and 
this approximation holds for the imaginary part of the 
LDOS as the surface is approached at the SPP frequency; 
the real part deviates significantly in this limit. However, 
when the incident frequency is detuned from the SPP 
(ui = 2.63 eV), the quasi-static approximation again be- 
comes valid. Similar results are found for the NLDOS, 
except when the inter-particle separation is greater than 
0.07 A and the quasi-static approximation again breaks 
down. This means that for resonance interactions one 
must be very careful about applying a quasi-static ap- 
proximation. 

It is also useful to examine the LDOS as a function 
of frequency for fixed particle position, as is shown in 
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FIG. 5. (Color online) Silver half-space system showing the coupled nano-particle forces, (a) Two dimensional graphs showing 
inter-particle forces as a function of height and particle separation; the particles are illuminated from above with a plane wave 
polarized in the y direction at the SPP frequency of the silver half-space and with the NP LSP tuned to be resonant with 
the SPP frequency. Colorbar indicates magnitude in log scale and arrows indicate directionality of the forces. One particle 
remains at xi = yi = and the second particle is moved in the y direction. The two particles are then both varied in the z 
direction such that they are always in the same plane. The arrows describe the force on the second particle which is not at the 
origin, (b) As in (a), but now the NP is far red detuned from the SPP frequency to a LSP resonance of u — 2.63 eV. (c) As in 
(a), but the incident field is at a frequency far detuned from the SPP frequency, ui = 2.63 eV. The NP LSP frequency of the 
particle is resonant with the SPP frequency, (d) As in (c), but with the NP ia also far red detuned to having a LSP resonance 
at lolsp = 2.63 eV. 



Figs. 4(a-b) for different NP detunings (see Fig. 2 and 
Eqn. 13). In both figures, the dashed lines correspond 
to p( Q \ the chain lines correspond to and the solid 
lines correspond to p^ 2 \ For simplicity we only focus on 
the imaginary part to examine the effects of the inter- 
actions. The particles both have their height fixed at 
h = 0.05 Xspp, and their separation is s = 0.05 Xspp- In 
Fig. 4(a) we consider a NP that is detuned by Acj = 1.1 
eV, and in Fig. 4(b) the NP is on resonance with the 
SPP. In Fig. 4(a), the SPP is visible at u = 3.73 eV in 
the LDOS when the particles are far detuned from this 
resonance, but the particles have a negligible effect on 
the SPP. Additionally, we see the resonances of the par- 
ticles interacting and produce a doublet feature caused 



by photon exchange effects. As the LSP resonances are 
moved towards the SPP resonance, the high frequency 
coupled LSP peak merges into the SPP resonance and 
acts to broaden it as well as detune it. The NLDOS be- 
havior (not shown) mirrors the effects seen here where 
the NPs strongly renormalize the NLDOS in the regime 
of the NP LSP regardless of where the LSP is with re- 
spect to the SPP. Similar effects for the LDOS and the 
NLDOS are seen in cavity- QED systems where the non- 
perturbative coupling between atoms or quantum dots 
causes additional photon exchange oscillations on top of 
the vacuum Rabi oscillations ' ' (the latter occur in sys- 
tems with suitably small dissipation). 

With the LDOS and NLDOS calculations acting as 
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a guide, we can now examine the light-induced force 
calculations for the geometry shown in Fig. 1 and de- 
scribed above. Four excitation regimes of interest shown 
in Figs. 5(a)-(d), corresponding to the regions highlighted 
in Figs. 4(a)-(b). We plot in log scale the intensity and 
use arrows to indicate the direction of the force. In 
Fig. 5(a) we illuminate at the SPP frequency and tune 
the LSP resonance of the particle to be at the same value. 
The particle separations and heights are both varied up 
to 0.25 A (= 83 nm). For an inter-particle separation 
greater than 0.15 A, the particles cannot feel each other 
except when the magnitude of the force is much less than 
1 pN, which happens at a height of w 0.125A and is due 
to the single particle interaction with the surface. The 
particles would thus be pushed away and then trapped in 
stationary positions w 0.125A above the surface and at 
a separation of 0.185 A. Interestingly, as the particles get 
closer to each other their interaction can still be negligi- 
ble compared to the particle-surface interaction if their 
height is smaller than their inter-particle separation; this 
is caused by quenching which reduces the transfer of ra- 
diation between the two NPs. If the inter-particle sepa- 
ration is sufficiently close, and greater than their height 
above the surface, then the particles strongly optically 
couple to each other and to the half-space - as can be 
seen by the fact that force still varies as the height of 
the particle varies. The vector force topology seen in 
this graph manifests itself through the coupling between 
the systems constituents as their separation varies. This 
dynamic coupling is very similar to the self-induced back- 
action demonstrated by Juan et al. 12. 

In Fig. 5(b), we consider a similar excitation scenario 
as in Fig. 5(a), except that the LSP of the NP is far red 
detuned, with ujlsp = 2.63 cV. In this case, we notice 
that the silver half-space dominates the response and the 
particles are continually drawn to the surface unless s < 
0.15 A. When s < 0.15 A the particles can couple to each 
other and are drawn together, however the effects of the 
surface seem to be negligible for h > 0.1 A. 

In Fig. 5(c) we tune our illumination to u> = 2.63 eV, 
but keep our LSP resonant with the SPP; note that 
0.25 A =118 nm. For particle separations greater than 
0.08A, the effects of the surface dominate, and for sep- 
arations below 0.08 A, the inter-particle effects dominate 
- though they still sensitively depend upon height. Both 
Fig. 5(b) and (c) exhibit very weak inter-particle cou- 
pling and very weak particle-surface coupling, so that 
the perturbative expression for Eq. (3) would hold. This 
is highlighted by the fact that the magnitude of the par- 
ticle forces are much lower than when we illuminate on 
the LSP resonance. 

For our final force example in Fig. 5(d), we examine 
the case when the LSP and the illumination are both 
far detuned (ui = 2.63 eV, cf lospp = 3.73 eV) from 
the SPP resonance. We observe three useful coupling 
regimes: i) When the particles arc very close to the 
surface (< 0.05 A), and the inter-particle separation is 
greater than 0.1 A, we see that the surface completely 



dominates the forces and the particles are pulled towards 
the half-space, ii) When the particles are very close to 
each other (< 0.1 A), then the inter-particle interaction 
dominates but this is again mediated by the half-space 
as there is a height dependence. Hi) In the remaining 
region, we can see that both inter-particle coupling and 
surface-particle coupling is present where the half-space 
dominates when h < s and particle-particle coupling is 
more dominant when h > s. This trend does not con- 
tinue indefinitely as for h > 0.2 A we see the inter-particle 
forces become weaker for equivalent separations. The role 
of electromagnetic quenching is also minimal as we are 
so far from the "lossy" SPP resonance. 

It is worth mentioning again, that the use of the gradi- 
ent force instead of Eq. (12) would predict an entirely dif- 
ferent answer. Additional calculations (not shown) show 
that for the case of Fig. 5(a), the gradient force topog- 
raphy is completely different, with an additional node 
along a vertical line at s/A 0.16 and no variation of 
force direction above hj A = 0.1. Thus using the gradient 
force for such a strongly perturbed system is generally 
not valid. 



IV. NEGATIVE INDEX MATERIAL SLAB 
WAVEGUIDE 

We next examine a 280-nm NIM metamaterial slab 
which supports a negative index in the frequency re- 
gion u = 0.78 - 0.92 eV. The relevant NIM and NP 
response functions are shown in Fig. 6. The possible 
benefits of using metamaterials is the ability to tune the 
material properties by engineering the constituents of the 
unit cell. Negative index metamaterials can be produced 
with very low loss in the microwave regime, however 
scaling to the visible has proven to be quite a challenge 
as the materials become very lossy ' , though continued 
progress is being made with new designs ''"". We will use 
NIM parameters that are close to experimental state-of- 
the-art for communications wavelengths, yet still have 
a respectable figure-of-merit: FOM = |Re(n)|/Im(n). 
Reported figurc-of-merits arc FOM = 2.0 at 1.8 /im 59 
and FOM = 0.5 at 780 nm 1 '"; for our calculations, 
FOM w 1.0 at u> = 0.78 cV. The permittivity is given by 
Eq. (13) with e r = 1, uj pe = 2.03 cV and 7 = 8.3 meV. 
The permeability is given by the Lorentz model , 



with the magnetic plasma frequency ui pm — 0.69 eV and 
the atomic resonance frequency ujq = 0.78 eV. 

Detailed descriptions of the exact corresponding com- 
plex band structure and Purcell effect are given by Yao 
et al. for the same parameters as given above. Here, 
we briefly point out a few features of interest for this 
study. At uj , the dispersion curves of all of the leaky 
slow light modes (SLMs) of the system converge at a sin- 
gle frequency, luslm = ^o; thus particles near the slab 
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FIG. 6. (Color online) Particle and NIM slab response func- 
tions. For both the NIM and the NP permittivity we use 
the Drude model [Eq. (13)], and for the NIM we model the 
permeability with a Lorentzian [Eq. (14)]. (a) Bare polar- 
izability [Eq. (1)] of a nanoparticle with the LSP resonance 
tuned to the SLM of the 280 nm thick metamaterial slab, 
uslm = ulsp = 0.78 eV (solid lines) (u; pe = 2.22 eV), and 
with the LSP resonance tuned to be off resonant with the LSP 
of the 280 nm metamaterial slab, ujlsp = 0.89 eV (dashed 
lines) (w pe = 2.53 eV). In both cases 7 = 11 meV. Gray- 
light curves indicates real parts and red-dark curves indicates 
imaginary parts, (b) Refractive index of the slab, where the 
gray-light curve indicates real part and the red-dark curve in- 
dicates imaginary part. Parameters are e r = 1, to pe — 2.03 
eV, 7 = 8.3 meV, cj pm 0.69 eV and u = 0.78 eV. 



can couple to many different slow light modes at this 
frequency. The slow light frequency regime gives rise 
to an enhancement in the LDOS and correspondingly to 
an increased Purcell effect. An increase in the LDOS is 
also seen near the SPP modes of the metallic surfaces, 
but the effects of quenching in metallic systems reduces 
the amount of light that escapes to the far field and the 
typical propagation distances of SPPs are limited by the 
material loss. However, slow light modes could, in prin- 
ciple, propagate for much longer distances. Also, the 
typical scaling laws associated with the quasi-static ap- 
proximation are not reached, even very close to the slab-' 1 
(because of the strong magnetic resonance) . Thus for our 
study, we are never really in the quasi-static regime above 
a NIM slab and we must consider retardation effects. It 
is also worth noting that NIM slabs support both TE 
(transverse electric) and TM (transverse magnetic) SPPs, 
which is in contrast to metallic surfaces that only support 
TM SPPs. The SPP modes in NIMs will not be discussed 
here, as their general properties are similar to the SPP 
of metals and at higher frequencies (ugp P = 0.92 eV and 
u™ p = 1.43 eV) 20 . 

For the metamaterial system, we consider a particle 
with a scaled radius of 0.015As£a/ = 24 nm, and we tune 
our NP to be in the frequency regime of our peak LDOS 
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FIG. 7. (Color online) (a) Complex LDOS, 
and (b) NLDOS, pffl (r, r',w), for a 280 nm thick meta- 
material slab as a function of height at ui = ojslm- For 
the LDOS r = (0,0, h), and for the NLDOS r' = (0,0, h), 
r = (0, 0.05A, h). Gray-light curves represent the real part 
and red-dark curves represent the imaginary part. The dashed 



the chain line 



line has no particles [p y (r;co), p yy (r, r'; 
has one particle located at the position the LDOS/NLDOS 
is calculated [pi, 1 ' (r;w), p y X y (r, r';w) and ri = (0,0, h)}, and 
the solid line has two particles; one where the LDOS/NLDOS 
is calculated at the other at the same height as the first but 
separated by 0.05A in the y direction [p y 2 ^ (r; cu), p y 2 y (r, r';uj), 
ri = (0,0, h), and r2 = (0, 0.05A, h)]. Gray shaded region in- 
dicates region where particles would overlap the surface. 



associated with the SLMs (at wo, see Fig. 6); practically, 
such tuning may be achieved, e.g., by using nano-shcll 
structures 62 . Additionally, we reduce the NP damping 
rate to 7 = 11 meV to examine the SLM features which 
would otherwise be obscured. For the metamaterial slab, 
we expect large enhancements of LDOS at the slow-light 
modes frequency similar to Ref. [20], however it is not 
obvious what the inter-particle coupling effects will be, 
nor the role of inter-particle coupling from the waveg- 
uide modes. Similar to the metal half-space case [Figs. 3 
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(a-b)], we first examine the LDOS and the NLDOS in 
Figs. 7(a-b) as a function of height. In Fig. 7(a), the real 
(gray-light curve) and imaginary parts (red-dark curve) 
of the LDOS again diverge as the slab is approached 
[Eq. (A10)] when no particles are in the system, however 
there is a change in sign of the real part compared to the 
metallic case, indicating the Lamb shift would be a red 
shift instead of a blue shift. The imaginary part (Purcell 
factor) reaches a value of 100 at 0.02\slm = 31.7 nm 
compared to the metallic case which reaches a value of 
100 at 0.05 Xspp = 16.6 nm. Thus the metamaterial 
gives an equivalent enhancement at twice the distance in 
absolute units. Introducing a particle at the location of 
the LDOS [ri = r = (0, 0, h)] renormalizes both the real 
and imaginary parts of the LDOS at small h/X and re- 
moves some of the divergence behavior for small distances 
close to the slab - similar to the metallic case. We also 
see that the maximum of the real part is no longer lo- 
cated closest to the surface. The addition of the second 
particle [r 2 = (0, 0.05A, h)] increases the imaginary part 
to a constant for h > 0.08 A to almost exactly the same 
value as for the metallic case which is due to the fact that 
we are in the quasi-static limit for the homogenous inter- 
action between the particles and the slab no longer plays 
a role. The real part dips slightly below zero indicating 
that there can be either a blue or a red shift depending 
on the height of the particles and stays below zero for 
h > 0.025 A. 

For the NLDOS [Fig. 7(b)], the real part (gray-light 
curve) follows a very similar trend as in the metallic case, 
where the zero particle case (dashed curve) is reduced as 
the slab is approached but is finite. Including the first 
particle (chain) drastically decreases the real part and 
thus the virtual photon exchange and the second parti- 
cle (solid) further reduces it. At closest approach the 
real part is small but still greater in magnitude to the 
metallic case by a factor of 20. The imaginary part with 
no particles (red-dark dashed curve) qualitatively follows 
the metallic case however when a particle is included in 
the system (red-dark chain curve), instead of reducing 
the real photon transfer there is an increase. The in- 
clusion of the second particle (solid curve) reduces the 
effect again but we still are able to increase coupling be- 
tween the particles compared to the metallic case. This 
transfer can be further increased as it crucially depends 
on the metamaterial loss used in the effective permittiv- 
ity and permeability. Obtaining lower losses is possible 
by improving metamaterial fabrication techniques which 
would result in less lossy slow-light propagation modes. 

A comparison of the LDOS in terms of frequency for 
the NIM slab is shown in Fig. 8 for different NP dctun- 
ings. Again, the particles both have their height fixed 
at z = z' = 0.05 Xslm = 79 nm, and their separation is 
0.05 Xslm- In Fig. 8(a) we consider a NP that is blue 
detuned by Acu = 0.11 eV, and Fig. 4(b), the NP is on 
resonance with the slab SLMs. When the NP is detuned 
from the slow-light modes (non- resonant case), we see 
that there is still a large enhancement of the LDOS at 
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FIG. 8. (Color online) Im[p { y N) (r,w)] for a 280 nm thick 
metamaterial slab as a function of frequency at the posi- 
tion r = (0, 0, 0.05Aslm). The dashed line has no particles, 
Imfpy ' (r;oj)], the chain line has one particle, Imfp^ 1 ' (r;uj)], 
and the solid line has two particles, \m[py (r; uj)], at locations 
n = (0,0, 0.05 Xslm), and r 2 = (0, 0.05A S £m, 0.05Xslm )■ 
(a) The particle LSPs are blue detuned by Aoj = 0.11 eV 
compared to the slab SLM frequency (see Fig. 6). (b) The 
particle LSPs are at the slab SLM frequency. Arrows indicate 
the particular scenarios we are examining in force graphs: 'i' 
- Fig. 9(a), 'ii' - Fig. 9(b), 'iii' - Fig. 9(c) and 'iv' - Fig. 9(d). 



the NP resonance compared with the zero particle case, 
but this is essentially the homogeneous space coupling 
due to the particles and is only slightly altered by the 
presence of the slab. When the NP is tuned to the ujslai 
resonance there is a greater enhancement than in free 
space but the coupling is dominated by inter-particle in- 
teractions. If we include only one particle, it is evident 
that the inter-particle effects dominate the spectrum as 
the LDOS more closely follows the zero particle case. 

To illustrate the effect on light-induced forces, we 
again consider four different cases that are highlighted 
in Figs. 8(a-d). We illuminate with a plane wave at the 
slow-light resonance frequency of the metamaterial slab 
(wslm = 0.78 cV), first with the NP on-rcsonance [Fig. 9 
(a)], and then with the NP off-resonance [Fig. 9(b)]. 
We then illuminate off resonance (u> = 0.89 eV) tuning 
the nanoparticle to be on-resonance with the slow-light 
modes [Fig. 9(c)], and then to be off-resonant and at the 
same frequency of the illumination [Fig. 9(d)]. All figures 
show the log of the magnitude of the force in intensity 
scale and arrows indicate directionality. 

When the NP and the slow-light modes are both on- 
resonance with the incident radiation [Fig. 9(a)], we see a 
very similar situation as when the NP was on resonance 
with the SPP [Fig. 5(a)]. For s < 0.2 A there is a divi- 
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(c) Force lines on particle 2, at region-'iii' on Fig. 8 
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FIG. 9. (Color online) Metamaterial slab system showing the coupled nano-particle forces, (a) Two dimensional graphs showing 
inter-particle forces as a function of height and particle separation; the particles are illuminated from above with a plane wave 
polarized in the y direction at the SLM frequency of the slab {ujslm = 0.78 eV). The NP LSP is also tuned to be resonant 
with the SLM frequency. Here one particle remains at X\ = j/i = and the second particle is moved in the y direction. The 
two particles are then both varied in the z direction such that they are always in the same plane. The arrows describe the force 
on the second particle which is not at the origin, (b) As in (a), but now the NP is blue detuned from the SLM frequency to a 
LSP resonance of u — 0.89 eV. (c) As in (a), but the incident field is at detuned to uj — ojlsp = 0.89 eV. (d) As in (c), but 
with the NP ia also detuned to having a LSP resonance at u = 0.89 eV. 



sion along the line s ~ 2 h where below this line the slab 
dominates the forces, above this line the inter-particle in- 
teraction dominates the forces and around which we see a 
combination of the two. As the height gets above 0.1 A we 
see that this does not continue indefinitely and the inter- 
particle interaction becomes weaker and shorter ranged 
as the slab no longer enhances the coupling between the 
two. Particles that have a small initial separation will be 
pushed away from each other and to a height of w 0.08 A. 
Note the forces here are an order of magnitude greater 
than in the metallic case, which is mostly due to the po- 
larizability scaling with the particle size but these could 
in principle be further tuned by improving the loss in 
these structures. 

When the NP is tuned off-resonance from both ujslm 



and the incident frequency [Fig. 9(b)], the long range 
coupling is lost compared to Fig. 9(a). For s > 0.12A, 
the force is dominated by slab interactions and are con- 
tinually drawn to the surface of the slab. For s < 0.1 A, 
and h < 0.1 A we see that the particles and the slab are 
all interacting which results in the particles being pulled 
towards the slab and together however the interaction 
range is short. When h > 0.1 A the particles essentially 
only interact with each other and are mostly drawn to- 
gether. 

Figures 9(c-d) show light-induced force calculations 
with the incident radiation detuned from the slow-light 
mode frequency to oj — 0.89 eV. In Fig. 9(c) the NP LSP 
is tuned to be resonant with the SLMs and off-resonant 
with the radiation. For s > 0.1 A the particles are un- 
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FIG. 10. (Color online) As in Fig. 9(a) but with a low 
msterial loss, (7' = 7/ 10) 



affected by each other and are repelled from the slab, 
except when they are almost touching. For s < 0.1 A 
the slab dominates over the inter-particle interaction for 
h < 0.1 A which is in contrast to Fig. 9(b). For h > 0.1 A 
the slab interaction weakens and the particles begin to 
interact and repel each other. 

Next, we examine the case where both the radiation 
and the NP are detuned to 0.89 cV and away from the 
SLMs, shown in Fig. 9(d). We see that, similar to 
Fig. 5(d), there are three different interaction regions, 
i)h< 0.08 A and s > 0.1 A, ii) h < 0.08 A and s < 0.1 A, 
Hi) and h > 0.08 A. In the first region, the slab domi- 
nates the force by pulling the particle when it is almost 
touching and pushing the particle away when it is slightly 
higher h > 0.02 A until the point where the vertical force 
becomes negligible and the particles are attracted to each 
other at h « 0.08 A. In region ii), the particles are 
causing a dramatic renormalization of the Green func- 
tion which leads to very strong, position dependent par- 
ticle interactions which are pushing the particles away 
from the slab and each other until they get to h rj 0.1 A, 
s « 0.1 A. Finally, above h = 0.1 A the particle interac- 
tion dominates but is mediated by the height above the 
slab and causes the particles to essentially be repelled 
to a fixed separation of s « 0.08 A, as the separation in- 
creases the slab starts to draw the particles towards it 
again however the particles will still be pulled towards 
s w 0.08 A. 

To investigate the influence of metamaterial loss on the 
inter-particle forces, we show the same scenario as Fig. 
9 (a) where the LSP and radiation is resonant with the 
SLMs, but we now decrease the material loss by a factor 
of 10, thus 7 = 0.83 meV. We show the resulting force 
in Fig. 10, and sec a number of important differences. 
First, where there was once a fixed height at which the 
particles would be attracted to, this height now varies 
with inter-particle separation. Below this dividing line 
in the region where h < 0.08 A and s > 0.18 A there are 
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(a) NLDOS for a nominal material loss, 7 = 8.3 meV 
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(b) NLDOS for a lower material loss. 7 = 0.83 meV 

FIG. 11. (Color online) Complex NLDOS, p { V y ] (r,r',w), for 
a 280 nm thick and metamaterial slab with regular loss (a) or 
with low-loss (b) as a function of separation [r' = (0, 0, 0.05A), 
r = (0, s,0.05A)] at cu — luslm- Gray-light curves represent 
the real part and red-dark curves represent the imaginary 



the 



part. The dashed line has no particles [p y J (r, r'; 
chain line has one particle located at the position the NL- 
DOS is calculated [p y y (r,r';cj) and ri = (0, 0, 0.05A)], and 
the solid line has two particles; one where the NLDOS is cal- 
culated at the other at the same height as the first but sep- 
arated along the y direction [p y 2 ^ (r, r';w), ri = (0,0,0.05A), 
and V2 = (0, s, 0.05A)]. Gray shaded region indicates region 
where particles would overlap each other. 



still inter-particle forces where in the regular loss case 
these forces have since died away. Finally, in the regular 
loss case, as the particles are moved vertically, along the 
line s = 0.04 A we see that the particles are attracted to 
each other close the slab, h < 0.125 A, but are repulsive at 
higher distances. This is contrasted in the low loss case 
where there is a division at s = 0.045 A, below which 
the particles arc always attracted and above which the 
particles are always repelled. 

To further examine how the loss alters the long-range 
coupling effects, wc plot the NLDOS in Fig. 11 for 
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both regular and low loss metamatcrial slabs at a fixed 
height, h = 0.05 A, and vary the separation between 
r' = (0,0,0.05A) and r = (0,s,0.05A). In both cases we 
see that the real part (gray-light curve) diverges at low s 
when there are no particles in the system (dashed line), 
due to the homogeneous part of the Green function. The 
addition of particles once more renormalizes these val- 
ues. We also see that in the low loss case, the real part 
plateaus between s = 0.08 A and s = 0.0125 A, whereas 
in the nominal loss case this simply decays, so material 
loss has a large influence on the light-induced forces The 
imaginary part of the NLDOS (red-dark curves) varies 
slowly towards zero in the regular loss case, however we 
see that the imaginary part of the NLDOS in the low loss 
case oscillates around a value of -2, with a much larger 
amplitude. The increase in the NLDOS allows the par- 
ticles to couple much farther in Fig. 10 than in Fig. 9 
(a). Thus for decreasing material losses, the NPs can 
be coupled over longer distances where this coupling is 
mediated by the slow light waveguide modes. 

V. CONCLUSIONS 

We have introduced a theoretical formalism to com- 
pute the Green function response of small particles within 
the vicinity of multi-layered geometries. We have applied 
this theory to calculate the non-perturbative force inter- 
actions between two NPs in the vicinity of surface plas- 
mon polariton modes for a metal half-space geometry, 
and in the vicinity of a slow-light NIM waveguide. Both 
planar structures facilitate a large local density of states, 
and non-local photon interactions between the particles. 
We have found that both structures exhibit rich but sim- 
ilar force maps despite the different mechanisms for in- 
creasing the LDOS. When both particle and slab (metal 
and NIM) are on resonance with the incident illumina- 
tion, the particles will be pushed away from each other 



and pushed to a fixed height above the slab (Figs. 5(a) 
and 9(a)). Such an effect would aid in preventing the ag- 
gregation of NPs. When the particle and illumination are 
off resonance with the slab then the particles are pulled 
to a specific height and pulled towards each other up to a 
fixed distance (Figs. 5(d) and 9(d)) which would enable 
the creation of dimers. In all the other cases the most 
likely scenarios are the particles being pulled towards the 
slab or the particles being pushed away from the slab. 

For metallic surfaces, the material parameters are 
largely fixed, limiting some of the engineering available 
to such structures, however metamaterials in principle 
have the ability to have their intrinsic parameters tuned 
by changing the basic unit cell. Such tunability will al- 
low simplified geometries such as the planar structures 
to aid in the creation of long range optical forces for 
the trapping and localization of small NPs. The same 
structures also exhibit rich and fundamentally interest- 
ing QED phenomena offering applications for radiative 
decay engineering of embedded quantum light sources. 
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Appendix A: Planar Green function 

The Green function above a planar structure can be 
written in terms of its angular spectrum 42 which involves 
decomposing the wavevector k = k x ic + k y y + k z z into 
its various components and integrating over each contri- 
bution. The real-space homogeneous Green function 



Ghom (r, r') 



-<5(R) 



oo poc 



f homG 



^(x-x') + ik y (y-y') +l k z \z-z'\ dk J, 



oo J — oo 



SB ' ' S7r 2 C 2 6B J- 

where ks = uj^/sb/c is the wavevector in the background material, and the z-component is given by k 
(k% — k\ — ky) . The matrix fh om is given by 

1 



(Al) 



fhom — 



h 2 — h 2 —h h 



x rvy rCy -T^y^z 

-fk x k z ^kyk z kg k z 



(A2) 



where the upper sign is used when z > z' and the lower sign is used when z < z' . The scattered part of the Green 
function in a multilayer environment (no particles) can be written similarly in terms of s and p polarized contributions, 



G sca tt (r, r ) — 



scatt 



J,k x (x-x')+ik y (a-a') +ifez ( z_z ')dfc a ;dfe 



k 2 k x ky 
- k x ky k x 











(A3) 
(A4) 
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r p (k x ,k y ) 



h 2 



7.2 



*2) 



fcsfc 2 



kx kyk z 
knkl 



'kx (k'i + ky) —k y (fc 2 + fc 2 ) 



k 9 ' 

K) 

-(fc 2 + fc 2 ) /k z 



kx \k x 
k y (fc 2 



(A5) 



Here the matrices (or dyadics) f scatt and f scatt are given 
in terms of the reflection coefficients, W p , for s and p 
polarization above the multilayer. For the three-layer ge- 
ometry of a slab with height h considered here, the upper, 
background layer having £q = £\,^b = Mi, the middle 
layer having £2, fJ.2 and the lower layer having £3, [13 these 
reflection coefficients are 



s/p 



.s/p.s/p s/p 2i/3h 
s/P 1 l 12 ''21 '23 e 



r s/p s/p 2ij3h ' 



(A6) 



21 '23 



1 /2 

where /3 = ± {uj 2 £2^2 / c 2 — fc 2 — fc 2 ) is the z com- 
ponent of the wavevector in the middle layer when 

Rc (w 2 £2/Wc 2 ) > Rc(fc 2 +fc 2 ) 1/2 . The sign of 
/3 depends on whether or not the refractive index 
of the middle layer is positive (upper) or negative 

(lower) 03 - 64 . For Re (w 2 e 2 // 2 /c 2 ) < Rc (fc 2 + fc 2 ) 1/2 , (3 = 

1 /2 

z (fc 2 + fc 2 — uj 2 '£2^2 1 'c 2 ) for both positive and negative 
index materials. The single-layer reflection and transmis- 
sion coefficients are 



/1 j k{ 



f^jkiz ~\~ ^%kj z 

2fijki z 
IJ-jki Z + i-iikjz 



£jkj z £jkj z 
£jkiz "t" £%kj z 



'>.) 



£jkiz H~ &ikj z 



(A8) 



Solutions of Eqs. (Al) and (A3) for real dielectrics can 
be difficult due to poles close or along the path of inte- 
gration, however this can be solved by numerically inte- 
grating around the poles in the complex plane which lie 
in a known region 4 ' . For lossy NIMs the poles along the 
integration path are found to be in the lower part of the 
complex plane, whereas for lossy positive index materi- 
als the poles are located in the upper half of the complex 
plane 20 . The solution described by Paulus et. 0Z. 47 is 
more complicated for materials which are able to sup- 
port negative index modes and surface plasmons as the 
location of the poles in the complex plane are essentially 



given by the complex band structure of the material 1 ' 1 . 
The largest contributions to the Green function are no 
longer confined to the region where Re [u 2 £2^2/^) > 

1 /2 

Re (fc 2 + fc 2 ) and careful attention must be paid to 
the integrand. This is trivial for a small number of cal- 
culations but can be cumbersome when many locations 
are required. As an example, for a two particle force cal- 
culation at a single point, the above calculations required 
14 separate Green function calculations when employing 
the dipole approximation. 

When considering the Green function in the quasi- 
static approximation, the homogeneous space Green 
function' ''' is given by 



Ghom,QS (r,r') 



4tt£bR 3 V R 2 



3RR j 



(A9) 



where T is the unit dyadic (diagonal terms are unity and 
non-diagonal terms are zero). The total Green function 
above a half-space in the quasi-static approximation in- 
volves the direct contribution from Eq. (A9) as well as 
the scattered contribution from an image source located 
beneath the surface 55 , 



^2 — £\ 

G QS (r, r ) = G ho m,QS (r, r ) T ; G homiQS (r, r ) . 



£2 + £1 



(A10) 



Here the minus sign is for x/y directed dipoles and the 
plus sign is for z directed dipoles. The location of the 
image charge is given by r" which is related to r' via, 
x' = x" , y' = y" , and z' = —z" when the surface of the 
half-space is located at z = 0. The scattering part of the 
Green function is then given by, 



G scatt ,QS (r, r') = T 



£1 



G 



£2 + £1 



hom,QS 



(r,r"). (All) 
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